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CN , Abstract 



To process data obtained during interference experiments in high-energy physics, 
methods of spectral analysis are employed. Methods of spectral analysis, in which an 



autoregression model of experimental data is used, such as the maximum entropy tech- 
nique as well as Pisarenko and Prony's method, are described. To show the potentials of 
the methods, experimental and simulated hummed data are discussed as an example. 

^ An alternative for an increase in the energy of colliding particles, observed when 

studying the properties of a substance at short distances, is the observation and study 
of the fine effects of the bound states of atoms and elementary particles |l|]. Of special 
interest are the interference energies and phases of states sensitive to the behaviour of the 
potential at short distances and to external effects ||. Information on the bound states 
of systems is obtained in interference experiments [EJ. The data obtained during such 
an experiment are represented by oscillating relations the spectral parameters of which 
are assessed to cast light on the energy and phase characteristics of a system M. The 



reliability of the resultant information depends on the ability to distinguish periodicities 
in experimental data and to estimate their parameters. 

The mathematical theory and algorithmic apparatus, used to reveal latent period- 
^ icities, have an extensive background and various practical applications. Methods for 

identification of periodicities are subdivided into a) approximation methods that enable 
a person to approximate experimental data by a function which a priori agrees with the 
preset model and b) filtration methods that provide information on components and their 
parameters when a priori evidence is minimum. The spectral characteristics of a signal can 
be estimated more accurately by approximation methods than by filtration methods |J. 

The authors discuss a search for latent periodicities in differently generated data using 
methods based on an autoregression data model (ARDM) ||. 

ARDM-based methods have an intermediate position between the groups of data pro- 
cessing methods discussed. In these methods, approximation is represented by a postulate 
of local data connection and filtration by operations done to produce a spectrum and to 
identify components. Owing to their intermediate position, ARDM methods take the 
advantages of each of the groups described. Spectral parameters and resolution can be 
estimated more accurately by ARDM methods than by filtration methods on the basis of 
Fouriers discrete transformation (FDT) 

Unlike FDT algorithms in which the signal observed is split into harmonic components, 
some of ARDM methods (e.g. Prony's method) use decomposition into components in 
the form of damping oscillations. When processing interference experimental data, such 
a representation has been found to be more realistic than a conventional one because it 
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considers the lifetime of states. This tendency persisted in Wavelet analysis, where de- 
composition into soliton-like components made it possible to follow space-time variations 
in the spectral characteristics of signals ||. 

Oscillation phases can be estimated more accurately by Prony's method than by FDT 
method. This advantage is important, considering the informative capacity of phases in 
interference experiments. 

ARDM modelling was done in the present study using the spectral maximum entropy 



data technique (MENT), Pisarenko's method and Prony's method [10-12], simulated and 
experimental data being used as an example. 

Autoregression data model. ARDM of a random process postulates a local data 
connection: the current value of the function y(l), which approximates the process studied, 
is estimated from the previous data counts x(l — k) in the form of a linear combination: 

p 

y{l) = Y,<k)x{l-k) (1) 

k=l 

with depth p and parameters a(k). The parameters of a(k) contain all information on the 
spectral characteristics of the process [T3|. 

Deducting the value of the variable x(l) from both parts of equality (1), the error in 
the current value e(l) = x(l) is expressed as: 

p 

e(l) = ^a(k)x{l - k) (2) 

k=0 

with the parameter a(0) = — 1. Because the error in the current value is also represented 
as a linear combination of previous data counts, the parameters of a(k) can be regarded 
as the coefficients of a linear filter which predicts an error. 

An error in predicting the current error value can be minimized, for example, on the 
basis of the least squares principle. In the case of an independent signal and additive 
noise, the observed data supplied to the input of the error prediction filter are converted 
to white noise. It is an example of a so-called "bleaching" filter. 

Multiplying both parts of relation (2) by the complex conjugate value x*(l) and aver- 
aging for all observed data, we will get a system of linear equations (3) for the parameters 
of a(k): 

p 

R(m - k)a(k) = a 2 5(m) (3) 

k=0 

where a 2 is noise dispersion, S(m) is delta-function, R(m — k) is the autocorrelation 
function values of R(k), m — 0, 1, 2 . . .p. 



is corn- 



To determine the filter coefficients of a(k), Derbu-Levinson's algorithm [14 
monly used by consecutively estimating bleaching filter coefficients more accurately and 
calculating the power of the error predicted, beginning with a first-order filter and ending 
with p-th-order filter. The algorithm calculates the sequence (an, of), (a 2 i, a 2 2, c|)> ■ ■ ■ j 
(dpi, a p2 , • • • , a pp , a 2 ) using the formulas: 

a n = -R(1)/R(0); of = (1 - \a n \ 2 )R(0) 
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a kk 



k-1 



R(k) + J2 a k-iM k 



i=l 



14-1 (4) 



The rapidity and economy of Derbu-Levinson's algorithm are due to the fact that, 
unlike the algorithm used to solve systems of Gauss linear equation, it requires p 2 opera- 
tions instead of p 3 operations and uses the results of previous calculations. More efficient 
algorithms have been developed for the same purposes [JT5[ . 



Maximum entropy technique (MENT). As a spectral method, it was first pro- 
posed by Burg |IJ to maximize the entropy density functional of the spectral capacity of 
a process P{1) 

(n-l)/2 

- lnP(/)^max (5) 

J=-(n-l)/2 

when fulfilling the conditions of Winer-Hinchin's theorem for p + 1 known values of the 
autocorrelation function 

(JV-l)/2 

P(l)exp{-jAculkAt) = R(k), (6) 

;=_(jV-l)/2 

where Aoj and At are frequency and time intervals between spectrum and data counts, j 
is imaginary unit (j 2 = —1), k — 0, 1, 2, . . .p. 

For the Gaussian random process, the requirement of a maximum entropy functional 
is equivalent to the minimum functional of a predicted error. This variation problem with 
the bleaching filter coefficients of a(k) as Lagrange factors is solved using the formula: 



P(l) 



a 2 



v 

a ( k ) exp(-jAulkAt) 

k=l 



(7) 



which has a simple meaning: the power spectrum of a signal is found by dividing the 
output power of noise by the squared module of the spectral characteristics of the bleaching 
filter. The coefficients of a(k) are estimated in ARDM by solving a system of equations 
(3). 

MENT considers local data connection, estimates the power spectrum of noise, does 
not result in negative values in the power spectrum and has a better resolution as com- 
pared with FDT, but is inferior in resolution to Pisarenko's and Prony's methods. As 
MENT is equivalent to the least squares method, it gives an undisplaced estimate of 
the spectrum. The requirement of maximum entropy or minimum information at pre- 
set restrictions automatically excludes all alien frequency peaks from the estimate of the 
spectrum. For this reason, the spectral peaks in MEM that correspond to an anharmonic 
periodic histogram are expected to be weaker than in FDT ||T6f . 

MENT is efficient for express estimation of the frequency spectrum of random pro- 
cesses. It can also be employed to reveal latent periodicities. To show this, let us discuss 
the dependence of the ionization potential of an atom in basic state on the atomic number 
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of an element Jl7j (Fig. la) and the time dependence of the rate of water flow through 
the Solomensky Strait of the Petrozavodsk Bay in Lake Onega |1| (Fig. 2a). To produce 



frequency spectra, programme [19] was used in these examples. 

The frequency spectrum in Figure lb shows periods 8, 10, 18 and 32 characteristic of 

Mendeleyev's Table. The FDT spectrum of such an anharmonic distribution contains, in 

i — 1 

addition to the main periods T k i, the periods T ki = T k i_x equal to 4.0, 2.7 and 2 for 

i 

the most intensive and narrow peak 8. In this spectrum, poorly intensive peaks are far 
more numerous, and among the periods 6.0, 4.7, 4.1, 3.1, 2.8, 2.3, 2.2 and 2.1 only 4.1, 
2.8 and 2.1 can be considered to be close to anharmonic periods. 

In the second example, the pattern of the signal studied is fairly close to a combination 
of harmonics and a noise component. In the power spectrum of this signal (Fig. 2b), all 
medium- and high-frequency peaks were interpreted. These peaks were correlated in 
the order of increasing periods with the seiches oscillations of the water body connected 
through the strait with the Petrozavodsk Bay, those of the Petrozavodsk Bay and those 
of Great Lake Onega. Correlation with observed field data and simulation data gave 
consistent results Spectral analysis revealed a ca. 12 hour tidal period, unknown 



earlier for lakes, and other low-frequency periods that were hard to interpret. The spectral 
resolution of MEM proved to be quite sufficient for analysis of the results obtained. 

Pisarenko's and Prony's methods. Pisarenko's method, used to distinguish har- 
monic components from their combination with white noise, provides even higher spectral 
resolution than MENT. It became possible to increase spectral resolution in Pisarenko's 
method at least 1000 times as compared with that of MENT. In this method, the fre- 
quencies of components are calculated, whereas in FDT or MENT they are estimated 
visually from the positions of the peaks of the frequency spectrum. Frequency estimation 
accuracy and spectral resolution are preset by the accuracy of calculations. 

To show the algorithm used to estimate frequencies in Pisarenko's method, let us 
discuss the sinusoid x(l) = sin(Ql). The trigonometric identity 

sin(flZ) = 2 cos Q ■ sin(fi(Z - 1)) - sin(fi(Z - 2)) (8) 

provides a link between data: 

x(l) = 2 cos Q ■ x(l - 1) - x(l - 2) (9) 

The Fourier's transformation of relation (9) results in the relation: 

X(u)- (1-2COSQ-Z- 1 + z~ 2 ) =0, (10) 

where z is understood as exp(— jfl). The roots of the equation of second degree Z\ = z\ 
determine the frequencies Q and — Q in accordance with the expression 

Q = arctg(ImZj/ReZj). (11) 

In a general case, local data connection in Pisarenko's method is preset in the form: 

x(l) = -J2<k)y(l-k)+s(l), (12) 

k=l 
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and for the coefficients of a(k) a system of linear equations: 

2p 

R{m - k)a{k) = a 2 a(m) (13) 

k=0 

is obtained. 

The coefficients a(k) and noise dispersion a 2 are estimated by calculating the eigen 
numbers and eigen vectors of the matrix of autocorrelation functions R(m, k) ||20|| . Dis- 
persion cr 2 corresponds to the minimum eigen number of the matrix R(m,k). In this 
case, it is convenient to use the recursive expression Rc(l + 1) = c(7), which follows from 
(13), and to determine the vector c(l + 1) from the vector c(l) estimated from previous 
iteration. It is convenient to take the vector c(0) = [1, 1, 1, 1] for initial approximation. 
Several iterations make it possible to produce a vector close to c(oo). It is used to estimate 

. 2 c T Rc _ . c 
Amin = <y = and vector a is determined as a = . 

C C -Amin 

The structure of the matrix of the autocorrelation functions R(m, k) is such that if 
harmonics are present in the process studied, then the filter coefficients a(k) are valid and 
meet the requirement: a(k) = a(2p — k). Composing and solving an algebraic equation of 
2p degree with the real coefficients a(k) and a(0) = 1 

Y^a{k)z 2p - k = 0, (14) 

fc=0 

complex conjugate roots, equal in module to unit, are calculated. The roots determine 
the frequencies of harmonics in accordance with expression (11). 

The positive frequencies of harmonics fli, i = 1, 2, . . . ,p, are used to estimate spectral 
power frequencies P(i), solving a system of equations 

v 

p (i) cos{ttilAt) = R(l). (15) 

i=i 

The self-consistency of the results obtained is checked by estimating noise dispersion a 2 
and expressing it by expression: 

a 2 = R(0)-J2m- (16) 

k=l 

Using Pisarenko's method, one can estimate the dispersion of additive white noise. How- 
ever, obtaining non-negative estimates of the spectral density of harmonics is not guar- 
anteed. Like in MENT, the initial phase of harmonics is not determined in Pisarenko's 
method. 

This shortcoming is overcome in Prony's method in which Pisarenko's idea is general- 
ized for non-steady periodic processes, using damping oscillations as a basis for spectral 
decomposition. 

In Prony's method, the function y(l), which approximates the process, is expressed 
through the complex numbers b{ and Zi by the expression 

2p 

y(0 = J>'^> * = 0,1, 2,.. .,n-l, (17) 
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bi = Ai exp(j8i) 

Zi = exp(«j + jfliAt) 

where the parameters Ai, ctj, f2j and 9i are the amplitude, coefficient of damping, fre- 
quency and initial phase of the i-th component of the decomposition of reference data. 

Local connection between data is postulated like in Pisarenko's method, the coefficients 
a(k) are also calculated and used to compose polynomial (14). The difference lies in the 
fact that the roots of the polynomial are now represented by complex exponents Zi. 

Proceeding from the minimum least squares functional 

71—1 

5>(Z) - y(l)) 2 - min, (18) 

1=0 

the vector B = (bi, b 2 , b 3 , . . . , b 2p ), which is dependent on the relation 

B = (F*F)- l F*X, (19) 

— * 

is estimated. In the above relation, X is a data line, F is Van der Mond's matrix composed 
of the degrees of roots z i: F* is a transposed matrix. From this relation and the roots of 
polynomial Zi the parameters of decomposition components are estimated: 

Ai = \bi\; Qi = arctg(InL2j/ReZ() 

^ = (In \zi\)/At; 0< = arctg(Im&i/Re&0 (20) 
Finally, Prony's method gives a data estimate: 

2p 

s(t) =J2 A * exp(aif) • exp(j(^t + 00), (21) 
i=i 

its Fourier image: 

2p 

S(u) = AeMJOih s (22) 

^ (uj- uJi) - jai 

and an energy spectrum: 

E(u) = \S(u>)\ 2 . (23) 

It has already been noted that an advantage of Pronys method is a broader class of 
functions used as a basis of spectral decomposition. The formalism of the method allows 
the incorporation of exponentially falling and rising signals into this class. Therefore, 
methods of spectral analysis can be applied to monotonically varying and transitional 
random processes. The opportunity to calculate the amplitudes, coefficient of damping, 
frequency and especially phase of oscillations can be used to interpolate and extrapolate 
signals. In some cases, this helps forecast processes, which is a difficult problem. 

To simulate Pronys method, the authors wrote a programme in Fortran. Figure 3 
shows initial signals in the form of a single sinusoid (a) and a mixture of four different 
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harmonics with additive white noise (b) and the results of their reconstruction based on 
the parameters estimated by decomposition. For the sinusoid preset by the formula 

1.600 cos(2.021t + 1.0798) (24) 
with 15% additive noise the estimate 

1.594 cos(2.021t + 1.0930) (25) 

was calculated, and for a mixture of harmonics 

2.000 cos(0.5054t + 0.2670) + 2.000 cos(1.015t + 0.5422) + 
+ 1.000 cos(1.844t + 0.948) + 1.600 cos(2..465t + 1.3167) (26) 

with 6% noise the estimate 

1.932 cos(0.5054t + 0.2821) + 1.957 cos(1.015t + 0.5887) + 

+ 0.9216 cos(1.844t + 0.9716) + 1.636 cos(2.465t + 1.2530) (27) 

The above results are in good agreement with the reference signals. It has been shown 
in these and other computing experiments that Prony's method estimates frequencies most 
accurately and the phases of components and the coefficients of damping of components 
least accurately and that the spectral parameters of high-frequency components are more 
accurate than those of low-frequency components. 

In the ARDM methods, a priori information, preset in the form of local data connec- 
tion, is generalized, can be used to form a model of an oscillating signal and to determine 
its composition and the spectral parameters of components without any special a priori as- 
sumptions. The resultant parameters can be estimated more accurately within the model 
formed (in some cases, it is desirable), using other powerful approximation methods, e.g. 
the Monte Carlo method . 

Conclusions: 

1. Methods of spectral analysis are highly informative. They are employed in interfer- 
ence experiments to cast light on the characteristics of the bound state of systems, such 
as atoms, nuclei and particles. 

2. ARDM-based methods, used to estimate the spectral characteristics of signals, 
combine the advantages of approximation and filtration methods. 

3. The estimates obtained by MENT for the power spectra of real signals are infor- 
mative, and can be used to analyse harmonic signals and anharmonic distributions. The 
spectral parameters of oscillating signals, estimated using Prony's method, are highly ac- 
curate and allow to represent experimental data by the superposition of basic functions. 
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List of Figures 



Fig. 1. Dependence of the ionization potential of atoms in the main state on the 
atomic number of an element (a) and the MEM estimate of its power spectrum (b). The 
frequency peaks of the power spectrum are arranged in accordance with the periods in 
Mendeleyev's Table. 

Fig. 2. Time relation of the current velocity through the Solomensky Strait of the 
Petrozavodsk Bay, Lake Onega (a) and its MEM estimate of the power spectrum. Fre- 
quency peaks are corresponded by seisches oscillation periods of 22-60 min in Lake Log- 
mozero (the lake connected through the Solomensky Strait with the Petrozavodsk Bay), 
1 hr and 20 minutes to 2 hrs and 18 min in the Petrozavodsk Bay, 3 hrs and 8 min to 4 
hrs and 27 min in Lake Onega and 12 hrs and 24 min in the tidal period. 

Fig. 3. The results of the simulation of Prony's method on models: (a) - a sim- 
ple harmonic signal 1.600 cos(2.02H + 1.078) and its estimate 1.594 cos(2.021t + 1.0930), 
(b) - a mixture of harmonics 2.000 cos(90.5054t + 0.2670) + 2.000 cos(1.015f + 0.5422) + 
1.000 cos(1.844t + 0.9848) + 1.6 cos(2.465t + 1.3167) and its estimate 1.932 cos(0.5054t + 
0.2821) + 1.957 cos 91.015* + 0.5887) + 0.9216 cos(1.844t + 0.9716) + 1.936 cos(2.465t + 
1.2530). 
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